clc;

% % == 1. 导入数据 ==
% % 读取 CSV 文件，假设第一列为时间，第二列为信号值
% fileName = '30kA0.csv';  % 替换为你的文件名
% 
% % 检查是否有标题行（使用 'NumHeaderLines' 来跳过首行标题）
% % disp('读取数据...');
% data = readmatrix(fileName, 'NumHeaderLines', 1);  % 自动跳过第一行文本
% 
% % == 2. 提取时间和信号列 ==
% time = data(:, 1);       % 时间列（单位秒）
% signal = data(:, 4);     % 原始信号列
% 
% % == 3. 确定采样频率 Fs 和归一化截止频率 Fn ==
% dt = mean(diff(time));  % 时间间隔均值（假设采样均匀）
% Fs = 1 / dt;            % 采样频率 (Hz)
% 
% % 奈奎斯特频率
% Fn = Fs / 2;
% 
% % == 4. 滤波器设计：5阶低通 Butterworth 滤波器，截止频率设为有效频率的 0.1 倍 ==
% order = 5;              % 滤波器阶数
% Fc = 0.2 * Fn;          % 截止频率（例子中为有效频率的10%）
% 
% % 归一化截止频率
% Wn = Fc / Fn;
% 
% % 设计Butterworth数字滤波器 数字低通 Butterworth 滤波器系数
% [b, a] = butter(order, Wn);  % 直接设计数字滤波器
% 
% % 或者换成 FIR 滤波器（更平坦滤波特性，如移动平均类）:
% % N = 30;
% % cutoffFreq = 0.02;  % 归一化截至频率 [0, 1] 对应到 [0, Fs/2]
% % lowPassFIR = fir1(N, cutoffFreq);
% % filteredSignal = filter(lowPassFIR, 1, signal);
% 
% % 将滤波后的信号相位对齐（常用方法避免信号延迟），使用 filtfilt
% filteredSignal = filtfilt(b, a, signal);  % 零相位滤波，支持在非实时数据中前后过滤
% 
% % == 5. 画图显示原始与滤波后信号 ==
% figure;
% 
% % subplot(2,1,1);
% % plot(time, signal, 'b');
% % title('原始信号');
% % xlabel('时间 [s]');
% % ylabel('幅值');
% % grid on;
% 
% % subplot(2,1,2);
% % plot(time, filteredSignal, 'r');
% % title('低通 Butterworth 滤波后的信号');
% % xlabel('时间 [s]');
% % ylabel('幅值');
% % grid on;
% 
% % 或者以下方式：将原图和滤波放在同一个坐标轴对比
% plot(time, signal, 'b', time, filteredSignal, 'r--');
% legend('原始信号', '滤波后信号');
% title('滤波前后波形对比');
% xlabel('时间 [s]');
% ylabel('幅值');
% grid on;


% 1. 导入波形数据
data = readmatrix('30kA0.csv','NumHeaderLines', 1); % 安全方式
time = data(:,1); % 单列时间
originalSignal = data(:,4); % 电压等信号

% 2. 第一次滤波 —— Butterworth 滤除多数噪声
Fs = 1 / mean(diff(time));  % 自动从采样算频率
[b,a] = butter(5, 0.1 * Fs / (Fs/2)); % 归化解，低通设定
filteredSignal = filtfilt(b,a, originalSignal);

% 3. 第二次粗糙平滑 + 放大时间解析度
% ==============================
smoothedSignal = sgolayfilt(filteredSignal, 1, 11); % 多一轮低通 + 拟一阶导数
timeHighRes = (0:length(smoothedSignal)-1) / (length(time)) * (time(end)-time(1)) + time(1); % 精度时间轴

% 4. 找出信号符号转变点 (有效过零区间)
% ==================================
signChanges = diff(sign(smoothedSignal)) ~= 0; % 找时间转变
zeroIndices = find(signChanges); % 例如，这里就是 clean signal 的 transisition idx

% 5. 高精度插值重构原始信号时间与信号
% ==================================
fSignal = interp1(time, originalSignal, 'spline');
timeInterp = linspace(time(1), time(end), 50000);  % 超过5e4
signalInterp = fSignal(timeInterp);

% 再一次过零点，找到最精细点位置 tzMatchedList(i)
% 还可以在此时构造零点 id 稳健处理：
% 零点符号 vs梯度穿，加速计算。

zeroXingTrueList = false(size(signalInterp));
for i = 1:length(signalInterp)-1
    if (sign(signalInterp(i))~=sign(signalInterp(i+1))) && (~(signalInterp(i) == 0)) % 创建双向过零间隔 <!-- 空粒 -->
        zeroXingTrueList(i) = true; % only mark start point
    end
end
zcPts = find(zeroXingTrueList);    % 使用插值后时间构建准确零流逝一步？

% 6. 对这些插值零点前后 ±1 μs 窗口线性拟合求取整体节段斜率
% ===========================================================
dt_min = 1e-6;                           % 1 μs
slopes = []; ident_z = [];               % 列表存储

for k = 1:length(zcPts)
    nIdx_ = zcPts(k);
    tz = timeInterp(nIdx_);               % 当前零点时间秒数码

    % · 构建相对于tz 的 +/- 1us窗, 精确拾取当前窗口下的信号：
    windowMid = timeInterp >= (tz - d_min) & timeInterp <= (tz + dt_min);
    localTime = timeInterp(windowMid);
    localSignal = signalInterp(windowMid);

    if length(localSignal) < 3  % 判定窗口信号是否有缺
        continue;  % 跳过处理
    end

    % · 时间对比 以tz为中心 zoom 跨度拟一阶线
    pfr = polyfit(localTime - tz, localSignal, 1);
    slopes = [slopes; pfr(1)]; % 取一阶系数为斜率

    ident_z = [ident_z; tz];
end

slopesClean = slopes(abs(slopes) <= 1e8);
% THRESHOLD
zc_TimeGood = ident_z(abs(slopes) <= 1e8); % = 移除外界
